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FINITE  ELEMENT  MODELING  OF  SUSPENSION  FLOWS 


S.  R.  Subia 
M.  S.  Ingber 

Deptartment  of  Mechanical  Engineering 
University  of  New  Mexico 


In  recent  years  the  research  community  has  devoted  some  effort  to  the  modeling  of  par¬ 
ticle  suspension  flows  using  a  continuum  approach.  Phillips  et  al.  (1991)  have  presented 
a  phenomenological  constitutive  model  of  shear-induced  particle  migration  in  suspension 
flows.  Because  particle  suspensions  are  used  in  numerous  manufacturing  processes  of  di¬ 
rect  interest  to  SNL,  we  have  begun  work  in  trying  to  model  these  problems  numerically. 
An  existing  SNLA  two-dimensional  finite  element  computer  code,  NACHOSII,  has  been 
modified  to  incorporate  the  equations  which  govern  particle  migration  in  suspensions. 

In  the  past  there  has  been  but  one  successful  implementation  of  particle  migration 
models  in  the  context  of  finite  element  method.  This  prior  implementation  was  an  explicit 
formulation  for  a  pipe  flow  problem,  but  our  present  work  is  based  upon  an  implicit  for¬ 
mulation.  Having  implemented  the  Phillips  model  of  shear-induced  particle  migration,  we 
have  spent  some  time  during  the  past  few  months  trying  to  benchmark  the  code.  To  date 
we  have  looked  primary  at  simulations  involving  wide-gap  Couette  devices.  These  problems 
have  been  chosen  primarily  because  experimental  NMR  data  is  available  for  these  cases. 

Concentric  Couette  Simulations 

Two  benchmark  simulations  of  concentric  wide  gap  Couette  devices  were  performed. 
Details  concerning  the  problem  geometry,  material  and  flow  parameters,  and  analytic  so¬ 
lutions  can  be  found  in  Phillips  et  al.  (1991).  Because  it  is  not  possible  to  model  these 
problems  as  strictly  steady-state,  simulations  are  performed  as  transient  evolutions  from  an 
initial  configuration.  Figures  11  and  12  summarize  the  result  of  two  problems  studied  for 
initial  particle  volume  fractions  of  50%  and  55%.  The  results  of  the  numerical  simulation 
indicate  good  correlation  with  both  the  analytic  results  of  Phillips  and  appear  to  fit  the 
experimental  data  best  in  the  center  portion  of  the  flow.  The  material  and  flow  parameters 
for  the  models  were  taken  directly  from  Phillips  and  no  attempts  were  made  to  adjust  these 
values  to  provide  a  better  fit. 

Eccentric  Couette  Simulations 

A  benchmark  simulation  was  performed  on  a  wide  gap  couette  device  in  which  the  inner 
cylinder  is  offset  from  the  outer  cylinder  center.  Unlike  the  case  of  the  rotating  concentric 
cylinders,  no  analytic  solution  exists  for  the  steady-state  particle  distribution  in  the  problem 
of  rotating  eccentric  cylinders.  Details  concerning  the  problem  geometry,  material  and  flow 
parameters,  can  be  found  in  Phan-Tien  et  al.  (1995).  One  numerical  simulation  was 
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Figure  1:  Particle  distribution  in  concentric  Couette,  initial  volume  fraction  50% 

performed  for  an  eccentricity  ratio  of  1/3  and  an  initial  particle  concentration  of  50%. 
Computed  results  are  compared  with  the  results  of  experiments  of  experiments  performed 
at  SNLA  in  Figure  13.  The  best  agreement  between  computed  results  and  experiment 
occur  in  the  center  portion  of  the  wide  gap.  Again  no  attempt  has  been  made  to  adjust 
the  material  parameters  for  a  better  fit  with  experiment. 

Piston  Flow  Simulations 

One  of  the  common  ways  of  displacing  a  suspension  in  a  tube  is  with  a  piston.  As  a 
prelude  to  studying  some  more  complex  suspension  transport  problems,  we  have  modeled  a 
simple  problem  of  suspension  displacement  between  two  pistons.  Simulations  were  carried 
out  for  piston  travel  to  diameter  ratio  of  10.  Qualitatively,  the  computed  results  are  found 
to  bear  some  resemblance  with  the  experimental  NMR  images  near  the  pistons.  However, 
the  center  region  between  the  two  pistons  have  yet  to  be  imaged,  and  will  later  be  compared 
with  the  computed  results. 
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Figure  2:  Particle  distribution  in  concentric  Couette,  initial  volume  fraction  55% 


Position 


Figure  3:  Particle  distribution  in  concentric  Couette,  initial  volume  fraction  55% 
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Applications  of  the  HP  Workstation  to  Problems  Involving 
Large  Deformations  and  Material  Failure 

Prof.  H.L.  Schreyer 

Department  of  Mechanical  Engineering 
University  of  New  Mexico 

One  of  the  proposed  tasks  for  which  the  workstation  was  requested  involved 
modelling  the  evolution  of  microcracks  along  grain  boundaries  in  almnina.  The  initial 
phase  of  this  research  has  been  accomplished  and  the  results  of  the  research  are  contained 
in  the  paper  by  Wang  et  al.,  (1995).  The  ability  to  use  the  workstation  was  an  essential 
ingredient  in  conducting  the  numerical  part  of  the  research  in  which  it  was  predicted  that 
microcracks  can  develop  and,  with  continued  loading,  some  microcracks  close  up  while 
others  evolve  to  macrocracks.  Such  features  were  also  shown  in  the  experimental  phase 
of  the  research. 

In  joint  work  with  Professor  Sulsky  of  the  Applied  Mathematics  Group  at  the 
University  of  New  Mexico,  a  new  computational  algorithm  (material  point  method)  has 
been  developed  for  analyzing  the  large  deformations  of  solids  (Sulsky  et  al.,  1994,  1995, 
1996).  This  capability  has  opened  up  the  possibility  of  obtaining  numerical  solutions  on 
a  workstation  that  were  previously  only  possible  on  mainframe  computers.  Examples  of 
problems  being  considered  are  extrusion,  metal  rolling,  injection  molding  and  solidification, 
metal  cutting,  penetration,  perforation,  upsetting  and  high  velocity  impact.  Sandia  National 
Laboratories  provided  a  graphics  package  which  is  particularly  suitable  for  displaying  the 
results  from  the  material  point  method.  This  package  has  been  iustalled  on  the  workstation 
and  provides  the  University  with  a  rmique  capability  of  providing  state-of-the-art 
calculations  and  graphical  output. 

The  enhancement  in  research  capability  provided  by  the  workstation  has  provided 
several  opportunities  for  joint  research  with  national  laboratories  and  industry.  In 
particular,  a  woman  graduate  student  from  Los  Alamos  National  Laboratory  performed  a 
detailed  study  of  metal  cutting  for  a  Master's  thesis  (Sharp,  1995).  The  only  comparable 
analysis  has  been  conducted  by  a  Ph.D.  student  at  a  major  university,  an  indicator  that  the 
procedure  developed  at  the  University  of  New  Mexico  competitive  with  much  fewer 
resources.  The  same  algorithm  is  being  used  as  the  basis  for  studies  on  metal  rolling  and 
on  the  inflation  of  airbags  by  Ph.D.  students  from  Sandia  National  Laboratories.  Recently 
Los  Alamos  National  Laboratory  has  adopted  the  procedure  for  a  study  of  the  interaction 
of  viscoelastic  particles  under  high  pressures. 

Professors  Sulsky  and  Schreyer  recently  visited  personnel  at  the  research  laboratory 
of  the  Aluminum  Company  of  America  where  intense  interest  was  expressed  in  using  the 
procedure  for  studying  the  compaction  of  aluminmn  particles  under  high  temperature,  a 
problem  that  is  extremely  difficult  for  conventional  finite  element  and  finite  difference 
procedures.  It  was  also  pointed  out  by  ALCOA  that  there  is  a  great  need  for  a  procedxire 
such  as  this  in  the  design  of  extrusion  molds. 

At  the  University  of  New  Mexico,  the  material  point  method  is  being  currently  used 
on  the  workstation  to  investigate  penetration,  the  failure  of  a  certain  class  of  composites 
by  delamination,  and  to  study  the  microstructural  effect  of  fibers  on  the  enhancement  of 
ductility  in  ceramics. 
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Although  these  research  activities  mvolve  a  considerable  amount  of  basic  analysis, 
there  is  no  question  that  the  presence  of  the  workstation  has  been  of  immense  benefit  in 
providing  numerical  results  to  illustrate  the  applications  of  the  theory.  For  example,  one 
consequence  has  been  that  Sandia  National  Laboratories  has  loaned  us  another  workstation 
to  perform  microstructural  studies  on  failure  mechanisms  of  electrical  connections.  In  the 
future  we  see  the  need  for  large-scale  calculations  and  with  appropriate  support  hope  to 
avail  ourselves  of  the  Maui  High  Performance  Computing  Center  operated  by  the 
University  of  New  Mexico  for  the  US  Air  Force. 
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Optical  Propagation  Through  a  Round  Turbulent  Jet 

C.  Randall  Traman  and  Rick  I.  Zadoks^ 

Department  of  Mechanical  Engineering 
University  of  New  Mexico 


Introduction 

The  effect  of  large-scale  (or  coherent)  structure  on  optical  propagation  through  a  round 
turbulent  jet  has  been  studied  using  numerical  simulations  as  well  as  experimental  data.  A 
passive  scalar  in  the  computations  is  related  to  refractive-index  fluctuations,  while  temperature 
fluctuations  in  a  heated  jet  were  used  in  the  experiment.  Linear  stochastic  estimation  was 
employed  to  determine  the  influence  of  the  large-scale  structure  on  the  refractive-index  field 
and  thus  on  optical  phase  error  and  beam  deflection.  The  importance  of  sharp  gradients  at  the 
edges  of  the  large-scale  structure  was  demonstrated.  Linear  stochastic  estimation  based  on 
three  temperature  probes  has  successfully  reproduced  the  majority  of  the  experimental  beam 
deflection.  When  linear  stochastic  estimation  is  applied  to  an  instantaneous  realization  from 
the  large  eddy  simulation,  velocity  and  vorticity  at  two  locations  are  sufficient  to  predict  the 
phase  error  through  an  entire  ring  vortex. 

Approach 

The  effects  of  turbulent  structure  upon  optical  propagation  through  a  free  shear  flow 
are  being  studied  by  measuring  and  predicting  the  dynamics  of  a  passive  scalar  in  a  low- 
Reynolds-number  round  jet.  The  experiments  are  carried  out  at  the  Air  Force  Phillips 
Laboratory  (AFPL)  in  collaboration  with  the  AeroOptics  working  group  led  by  Dr.  Lenore 
McMackin. 

Linear  Stochastic  Estimation  (LSE)  has  been  used  to  analyze  coherent  turbulent 
structure  in  the  experimental  data  and  in  numerical  simulations  of  the  roxmd  jet,  for  both  the 
distorted  wave  front  and  field  variables.  The  capability  of  conditional  averages  formed  by 
LSE  to  reproduce  turbulent  structure  relevant  to  optical  phase  distortion  is  being  studied  for 
the  AFPL  data  and  a  low-Re  numerical  simulation  database.  The  coherent  structure  in  free 
shear  flows  is  characterized  by  sharp  interfaces  caused  by  entrainment  of  fluid  from  one 
stream  deep  into  the  other.  The  instantaneous  structure  of  the  refractive-index  field  (or 
passive  scalar)  is  produced  by  the  large-scale  turbulent  structure.  Thus  the  optical  distortion 
induced  in  an  optical  beam  propagating  through  the  flow  will  depend  directly  on  the  dynamics 
of  this  large-scale  structure.  Compensation  or  correction  of  the  distorted  wave  front  will  thus 
require  detailed  knowledge  of  the  flow  dynamics. 

The  criteria  upon  which  conditional  averages  are  based  are  being  examined  to 
determine  whether  LSE  based  on  limited  experimental  information  can  adequately  represent 
the  scalar  field  to  predict  the  optical  distortion  across  the  flow.  Such  data  will  include  a  few 
point  measurements  of  velocity  or  temperature  and  the  deflection  of  one  or  more  thin  beams 
propagating  across  the  shear  layer  which  are  provided  in  the  AFPL  experiments.  It  has  been 
confirmed  that  these  few  measurements  can  provide  enough  information  to  closely 
approximate  the  large-scale  features  of  the  optical  phase  distortion. 


^Now  at  Mechanical  Engineering  Department,  University  of  Texas,  El  Paso. 
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Progress 

We  have  developed  an  experimental  turbulent  jet  facility  at  the  Air  Force  Phillips 
Laboratory  in  which  flow  and  optical  parameters  can  be  measured  simultaneously  (Truman,  et 
al.  1994).  Acoustic  forcing  was  added  to  the  plenum  chamber  to  enhance  vortex  formation  in 
the  jet.  Measurements  of  beam  deflection  through  a  turbulent  jet  for  various  beam  diameters 
(or  apertures)  showed  that  maximum  deflection  occurs  for  the  smallest  beam  diameter  and 
verified  the  importance  of  sharp  gradients  at  the  edges  of  large-scale  turbulent  structures 
(Truman,  et  al.  1994).  We  have  also  implemented  Linear  Stochastic  Estimation  (LSE)  of 
beam  deflection  for  propagation  through  the  center  plane  of  the  experimental  jet  flow  based 
on  3  temperature  probes  within  the  jet  shear  layer  (Truman,  et  al.  1995).  Data  from  a  large 
eddy  simulation  of  a  turbulent  round  jet  which  includes  a  passive  scalar  was  acquired.  LSE 
of  phase  error  and  beam  deflection  for  propagation  through  the  center  plane  of  the  jet  based 
on  2  reference  points  within  the  jet  structure  where  velocity  and  vorticity  components  are 
known  was  implemented  (Barsun,  et  al.  1996). 

Linear  Stochastic  Estimate  (LSEl.  Linear  Stochastic  Estimation  (LSE)  is  a  simple  yet 
powerful  technique  for  approximating  conditional  averages  of  correlations  between  turbulent 
quantities  (Adrian  et  al.  1989).  The  utility  of  conditional  averages  is  that  turbulent  events 
(e.g.,  intermittent  bursts)  can  be  accurately  identified  with  limited  flow  field  information,  such 
as  velocity  fluctuations  at  a  few  predetermined  points  in  the  flow.  The  LSE  is  formed  from 
two-point  spatial  (or  temporal)  correlations  between  the  quantities  of  interest,  such  as  velocity 
and  scalar  fluctuations  or  even  optical  phase  error.  The  LSE  of  a  turbulent  shear  flow  gives  a 
remarkably  good  reproduction  of  the  large-scale  structure  (e.g.,  Glauser  et  al.,  1993). 
According  to  Adrian  et  al.  (1989),  the  "imconditional  spatial  correlations  have  imbedded 
withm  them  the  structures  of  the  conditional  fields." 

LSE  of  beam  deflection.  Three  cold- wire  temperature  probes  were  placed  in  the  jet 
shear  layer  immediately  downstream  of  a  thin  laser  beam  propagating  through  the  flow.  After 
passing  through  the  heated  jet,  the  beam  is  focused  on  a  lateral  effects  detector  (LED)  to 
measure  its  deflection.  The  probes  as  well  as  the  beam  lie  on  the  center  plane  of  the  flow; 
only  streamwise  deflections  are  considered  here.  Probes  are  placed  on  either  side  of  the  jet  to 
include  the  influence  that  each  side  of  the  shear  layer  has  on  the  beam  deflection. 

Spectra  from  the  probes  and  the  beam  deflection  were  used  to  select  a  forcing 
frequency  which  excited  the  shear  layer  mode  of  the  nozzle  exit  boundary  layer.  Acoustic 
forcing  provided  by  a  small  audio  speaker  in  the  plenum  was  sufficient  to  enhance  vortex 
formation  at  this  fundamental  frequency. 

Temporal  correlations  between  the  beam  deflection  and  each  temperature  signal 
illustrate  the  differences  in  the  structure  of  the  temperature  field  detected  at  each  probe 
location.  Beam  deflection  in  the  streamwise  direction  is  directly  related  to  streamwise 
gradients  in  the  temperature,  which  vary  as  vortices  convect  through  the  measurement 
locations.  The  LSE  was  created  using  the  correlation  between  beam  deflection  and  each  of 
the  three  temperature  signals  with  appropriate  time  delays.  The  characteristic  behavior  of  the 
beam  deflection  was  found  to  be  predicted  by  the  LSE  very  well.  We  are  studying  the 
selection  of  different  time  delays  in  the  correlation  as  well  as  different  probe  locations  at 
several  streamwise  locations. 

Large  Eddy  Simulation.  We  have  examined  the  Large  Eddy  Simulation  (LES)  of  the 
low-Reynolds-number  turbulent  jet  (Chen  et  al.  1993)  to  determine  the  relationship  between 
the  scalar  fluctuations  and  the  velocity  and  vorticity  fields.  While  the  scalar  fluctuations  can 
be  linked  to  streamwise  velocity  fluctuations,  the  scalar  distribution  produced  by  the  rollup  of 
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the  ring  vortex  is  complex.  Some  geometric  characteristics  of  the  scalar  distribution  must  be 
linked  to  azimuthal  vorticity  or  radial  velocity.  Two-point  correlations  between  the  scalar  and 
these  three  variables  were  computed  within  the  first  highly  developed  ring  vortex.  Each 
correlation  was  computed  from  an  ensemble  of  values  from  around  the  jet  at  fixed  azimuthal 
displacement  A0  for  each  (r,Ar,z,Az)  since  the  flow  is  inhomogeneous  in  r  and  z. 

The  scalar  field  in  each  plane,  0  =constant,  was  computed  by  a  Linear  Stochastic  Estimate 
(LSE)  with  two  fixed  reference  points  in  the  plane  0=0.  Streamwise  velocity,  radial 
velocity  and  azimuthal  vorticity  at  these  reference  points  were  the  input  to  the  LSE.  The 
estimates  of  the  scalar  field  were  seen  to  be  quite  good  over  the  entire  spatial  extent  of  this 
vortex  (Barsun,  et  al.  1996). 
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